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Two coupled BECs with a large population imbalance exhibit macroscopic quantum self-trapping 
(MQST) if the ratio of interaction energy to tunneling energy is above a critical value. Here we 
investigate effect of quantum fluctuations on MQST. In particular, we analyze the dynamics of a 
system of two elongated Bose gases prepared with a large population imbalance, where due to the 
quasi one dimensional character of the gas we no longer have true long range order, and the effect of 
quantum fluctuations is much more important. We show that MQST is possible in this system, but 
even when it is achieved it is not always dynamically stable. Using this instability one can construct 
states with sharply peaked momentum distributions around characteristic momenta dependent on 
system parameters. Our other finding is the nonmonotonic oscillating dependence of the decay rate 
of the MQST on the length of the wires. We also address the interesting question of thermalization 
in this system and show that it occurs only in sufficiently long wires. 



I. INTRODUCTION 

The Josephson effect [H 13, S and macroscopic quan- 
tum self trapping (MQST) [J,|5| have been studied exten- 
sively both theoretically and experimentally in the con- 
text of two coupled BECs. MQST has been observed ex- 
perimentally in a single bosonic Josephson junction Q. 
In these cases, one considers a two well system prepared 
with a large initial population imbalance and investigates 
the dynamics. The occurrence of MQST can be under- 
stood by mapping the system into an equivalent classical 
nonrigid pendulum system, with the pendulum length 
depending on angular momentum Q. Angular momen- 
tum in the pendulum maps to population difference in the 
coupled BEC, and likewise angular displacement maps to 
the relative phase. Then the Josephson effect is analo- 
gous to small oscillations of the pendulum around equilib- 
rium (where the system has low energy), and MQST cor- 
responds to the pendulum making full revolutions (high 
energy). In the high energy pendulum case, the average 
angular momentum is nonzero and the angle monotoni- 
cally increases in time. This corresponds to nonzero av- 
erage population difference and a running phase in the 
BEC case, which defines MQST. 

MQST can be understood using the following qualita- 
tive considerations. Consider the case where initially all 
atoms occupy one well. The transition to MQST is a non- 
linear effect governed by the parameter A = J//i, where /i 
is the chemical potential in the full well (which is approx- 
imately equal to its interaction energy per atom), and J 
is the tunneling energy between the two wells. Self trap- 
ping occurs when A is smaller than some critical value of 
the order of unity. The transition to MQST is understood 
in terms of a simple energy argument. When the param- 
eter A is small, corresponding to large interaction energy, 
the initial state with a large population imbalance has a 
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very high energy. We might expect the system to dynam- 
ically equilibrate and the imbalance to disappear. How- 
ever, there is a kinematic barrier for this to occur since 
the interaction energy has to be transferred to the kinetic 
energy of atoms. If the condensate depletion is small and 
the dynamics is effectively constrained to the condensate 
modes (this is the case e.g. for weakly interacting 3D con- 
densates) then the only source of the kinetic energy is the 
tunneling term. For small A the gain in hopping energy 
cannot make up for the loss in interaction energy thus 
there is a phase space restriction to dynamics. The best 
the system can do is a partial tunneling of atoms, and 
the system is constrained to remain imbalanccd. This is 
analogous to giving the pendulum a strong enough initial 
kick to make full revolutions, where the angular momen- 
tum remains nonzero at all times. Here we have a large 
average angular momentum, where deviations from this 
average (due to the presence of the external gravitational 
field) oscillate in time with a high frequency and small 
amplitude. It is these angular momentum deviations in 
the pendulum case that are analogous to the high fre- 
quency small amplitude tunneling of atoms in the bosonic 
Josephson junction in MQST Finally, let us briefly 

mention that we can also understand MQST through an 
analogy with the emergence of second order atom tun- 
neling in double well systems due to superexchange (see 
Rcfs. @,0| and the endnote [11]) 

The explanation above of MQST heavily relies on the 
assumption that the condensate remains coherent at all 
times. This assumption is essential since the interac- 
tion energy can be also transferred to the kinetic energy 
of the non-condensate modes within each well. As we 
mentioned this process is kinematically suppressed if we 
are dealing with three-dimensional condensates where the 
depletion is small and quantum fluctuations required to 
excite such modes are negligible. In one-dimensional sys- 
tems, on the contrary, the situation is quite different. 
Quantum fluctuations are very important. In equilibrium 
they result in destruction of the condensate in the ther- 
modynamic limit even at zero temperature @. One can 
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expect that these fluctuations will also affect dynamics. 
This is indeed the case in a different setup, where it was 
shown that quantum fluctuations were responsible for de- 
struction of the effective spin echo in one-dimensional 
condensates Q . As we will explain in this paper the non- 
condensate modes crucially affect the dynamics expected 
from the simplified self-trapping picture. 

In this paper we will consider the two wire system, 
where all atoms are initially loaded into one wire and 
then the wires are suddenly coupled. Such a system was 
recently realized experimental ly I l0l| and analyzed theo- 
retically in a different regime We assume that ini- 
tially the atoms are in the ground state of a single full 
wire, and then the system is quenched by suddenly cou- 
pling the wires by the tunneling term. In this work we 
will focus on the self-trapping regime where the tunnel- 
ing is much smaller than the interaction energy per atom 
in the full wire. This initial state then describes a highly 
nonequilibrium situation and we can expect interesting 
dynamics. We then ask, is it possible to have MQST, 
and is it stable? If not, what is the mechanism of MQST 
decay? Does the system eventually equilibrate, and if so, 
at what time scales? 

To tackle the problem we will use the idea that in 
the self-trapping regime we have clear separation of time 
scales in the dynamical process. The high frequency 
small amplitude oscillations expected from zero dimen- 
sional analogy with Josephson junction [1,01 j will then 
serve as a source term for slow dynamical processes occur- 
ring at much longer time scales. The whole dynamics is 
then analogous to that in an externally driven system by 
a periodic modulation with frequency much larger than 
the natural transition energy in the system. Then one 
can effectively integrate out high frequency modes by av- 
eraging dynamics over short time scales and get coarse 
grained effective description of the slow degrees of free- 
dom (see e.g. Ref. [iJl for illustrations of this idea in 
the context of classical mechanics). This analogy allows 
one to construct an effective "low energy" description of 
the dynamics in our system describing slow degrees of 
freedom. The analogy with externally driven systems is 
not perfect, however, since in our situation there is a 
feedback from the slow modes affecting the fast coherent 
dynamics. As we will show this process inhibits original 
fast oscillations and leads to decoherence and eventual 
thermalization in the system. We have to rely on numer- 
ical simulations to describe this process. Nevertheless at 
intermediate times the picture where one treats high fre- 
quency degrees of freedom separately providing external 
source to low frequency fields gives very accurate quanti- 
tative description of the dynamics and allows us to derive 
analytical expressions describing key instabilities in the 
system destroying self-trapping. 

Another important result of the present analysis is the 
evidence of long-time thermalization in the system. We 
note that the dilute interacting ID Bose gas is well de- 
scribed by a delta function interaction, which is inte- 
grable in both classical and quantum descriptions p^ . 



and therefore does not thermalize in the usual sense [ij] . 
However, the coupling between the two ID wires breaks 
integrability and we might expect the system to thermal- 
ize. The issue of thermalization in weakly nonintegrable 
systems is far from clear (see e.g. Ref. [lal)- In a classi- 
cal system, the KAM theorem states that adding a weak 
perturbation to an intcgrable systern precludes the sys- 
tem from thcrmalizing (see Refs. [H [13; for the original 
KAM papers see Refs. [lllilllSl). It is very little known 
what happens in quantum weakly integrable systems. In 
this work we address this issue numerically for our partic- 
ular problem and give evidence that like in the classical 
systems thermalization may or may not occur in the an- 
alyzed setup depending on the choice of parameters. We 
note that the similar conclusions of thermalization due 
to breaking integrability were reached in Ref. [21| about 
thermalization in a single wire due to virtual excitations 
of higher radial modes. This paper concluded that ther- 
malization should occur for arbitrarily weak coupling to 
virtual excitations. However, experiments show that non- 
thermalizable dynamics in quasi one-dimensional systems 
of bosons can be very robust p^ . so the whole picture 
of thermalization in nearly integrable quantum systems 
remains quite controversial. 

Microscopically the system we are considering is de- 
scribed by the Hamiltonian (here and throughout the 
paper we set h = 1) 

-J dx (^-l^-a + ^^^i) , (1) 

where i = 1,2 denotes one of the two wires, g denotes 
the interaction parameter, m is the mass of the atoms, L 
is the system size, and J is the tunneling matrix element 
between the two wires. All the parameters are positive, 
and for simplicity we will use periodic boundary condi- 
tions. The operator '^i{x) annihilates an atom located 
at X on wire i. Wc assume that initially one of the wires 
is occupied by atoms with the density po while the other 
wire is initially unoccupied. At t = we couple the wires 
with finite, but small J so that we mimic the conditions 
that give rise to the MQST in the two well case. The 
initial chemical potential in the full wire /x = gpo is sim- 
ply the interaction energy per atom. Like in the two 
well case, the parameter A = J/gpo = J/p-, which is the 
ratio of interaction and tunneling energies, controls the 
transition between MQST and no MQST. There is also 
a natural length scale in this system equal to the healing 
length of the full wire ^ = 1/ y'mgpQ ~ l/^Jmp. Simi- 
larly the chemical potential p = gpo defines the (inverse) 
natural time scale in the system. In the natural units 
where r = gpot and x is given in terms of ^, the speed of 
sound in the full wire is = 1. 

The static and dynamical properties of the system are 
completely governed by the following three parameters: 
£ = L/£^, K = TrpoC: ^- Iii the weakly interact- 
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ing regime the parameter K ^ I is equal to the effec- 
tive Luttinger liquid parameter for the full wire Q. In 
equilibrium K sets the scale for quantum fluctuations 
in the system. As interactions increase the expression 
K = 7rpo$ for the Luttinger liquid parameter is no longer 
valid and should be substituted by another function of 
Po^i which can be extracted from the Bethe ansatz so- 
lution [23j . In the regime of infinite interactions (Tonks 
gas) the Luttinger liquid parameter saturates at K = 1. 
In this work we will focus on the regime of weak quantum 
fluctuations {K/tt > 10) corresponding to the semiclas- 
sical limit. As we will show in this (weakly interact- 
ing) regime the dynamics is very nontrivial. Our goal 
is then to treat quantum fluctuations perturbatively and 
analyze how they destroy the purely classical MQST. In 
this limit one can study the Hamiltonian ( [T|) both an- 
alytically and numerically usin g th e semiclassical trun- 
cated Wigner method (TWA) M, M, [13, [H, HI- 
This method gives the leading quantum corrections to 
the classical (Gross-Pitaevskii) dynamics accounting for 
the zero point fluctuations of the bosonic creation and 
annihilation fields in the initial state. The inverse of the 
Luttinger liquid parameter plays the role of the effective 
Planck's constant. It can be shown that the higher order 
quantum corrections to TWA appearing in the form of 
quantum jumps are down by a factor of [2^ . In the 

regime of interest such corrections are negligible so TWA 
is expected to give accurate quantitative description of 
quantum dynamics in our system. 

We note that in an infinite system for a ID gas there 
is no true long range order and no true condensate. At 
best we have a quasicondensate at zero temperature, so 
we have a power law decay of correlation functions or 
the off-diagonal elements of the single-particle density 
matrix: (l'^(x)l'(O)) l/\x\^/'-^^\ However, at finite 
size systems one still has a condensate fraction since a 
single mode of the system is macroscopically occupied 
compared to all other modes. Unlike in the 3D case, the 
condensate fraction strongly depends on the system size. 
But if we confine ourselves to small enough £ and large 
enough K so that 

1 /n'/'"^ 

is satisfied then we have a small condensate depletion and 
thus a highly coherent gas. Physically Eq. ^ gives us a 
rough estimate of the condensate fraction, and if Eq. ([2]) 
is satisfied we can confidently say we have a condensate. 

Provided we satisfy the condition ( [2]), and we have 
small enough A (as in the bosonic Josephson junction 
case), we show that it is possible to achieve MQST even 
in the ID case. By using an effective model descrip- 
tion, we explicitly find the conditions that the system 
must satisfy in order to obtain MQST, and we numeri- 
cally confirm these results. We also find that MQST is 
generally unstable, and the system eventually (on a time 
scale on the order of 1/J) decays from MQST through a 
dynamical instability. Nevertheless, the system exhibits 



characteristic MQST behavior (i.e. high frequency small 
amplitude tunneling of atoms between the two wires) on 
time scales less than 0(1/ J). More interestingly, with a 
judicious choice of the system parameters £, A, and K, 
we can use the decay mechanism to resonantly populate 
characteristic momentum modes Qc in the initially empty 
wire and thus get a very sharp momentum distribution 
at intermediate times. Depending on the choice of pa- 
rameters, the system may or may not thermalize. We 
will show what is necessary in order to achieve this by 
numerically solving for the dynamics for the Hamilto- 
nian (see Eq. ([T])). Through our effective model, we can 
also attain a quantitative understanding of the underly- 
ing mechanism leading to the resonant transfer of atoms 
to characteristic modes and subsequent decay of MQST. 

The paper is organized as follows. We first discuss the 
approximate effective description of the system dynamics 
in Section |TT1 This description gives us the mechanism 
for the dynamical instability and a quantitative under- 
standing of how to achieve the sharp momentum distri- 
bution at intermediate times. In Section IIIII we show 
specific results from the numerical solution of the micro- 
scopic Hamiltonian ([T]), detailing what is necessary for 
MQST, the question of MQST stability and the forma- 
tion of the intermediate sharp momentum distribution 
state, and address the question of thermalization of the 
system at later times. Finally in Section IIVI we discuss 
and give an interpretation of our results. 



II. EFFECTIVE MODEL 

Given the microscopic Hamiltonian in Eq. ([1]), consider 
initial conditions where the two wires are initially decou- 
pled and all atoms occupy the first wire {i = 1). Also 
assume that the first wire is prepared in such a way that 
it is close to the ground state (T ~ 0), and take the Lut- 
tinger K and system size L such that the depletion from 
the zero momentum mode is small, i.e. assume that the 
condition ^ is fulfilled. Then in general for short times 
(< ^ 1/J) after the quench, the dynamics will be domi- 
nated by zero momentum modes as it is in the two well 
systems where MQST was previously studied. 

As explained in the introduction, this occurs because of 
the large separation of scales present in our current sys- 
tem if we choose system parameters such that we are in 
the MQST regime. In particular, we anticipate that the 
zero momentum (condensate) mode will undergo small 
amplitude high frequency oscillations between the two 
wires. We are then justified in treating the high frequency 
dynamics separately from long time dynamics governing 
non condensed modes. Thus in the leading order we can 
treat ID systems as zero dimensional system and analyze 
condensate dynamics separately. Because of the high oc- 
cupation of the zero momentum mode in the full wire, we 
have to use a full nonlinear treatment of the dynamics. 
This can be easily done in this case, however, as we only 
have to treat a single mode. This mode can be treated 
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classically (using Gross-Pitaevskii approach). In the next 
order of approximation we will treat this fast condensate 
dynamics as an external source for the non-condensate 
modes. We can anticipate that the momentum modes, 
which are resonant with the fast condensate oscillations 
should be excited. As we will show such process can be 
viewed as enhancement of the zero point fluctuations of 
high-momentum modes. Therefore we have to treat such 
modes quantum-mechanically. At short times while the 
occupation of the non-condensate modes is small we can 
linearize the resulting equations of motion (with time- 
dependent parameters) and obtain analytic results. At 
longer times excitations of these modes should provide 
damping to the condensate oscillations and eventual de- 
struction of the condensate dynamics. This is a nonlinear 
effect happening at longer time scales, associated with 
the tunneling coupling 1/J (as opposed to the period of 
fast oscillations l//x). In the MQST regime we have a 
large separation of scales because the ratio A = J/jU is 
small. So our strategy to address the problem is justi- 
fied. 



A. Effective Hamiltonian for zero momentum 
modes 

Our starting point in the analysis of this system is to 
take the Hamiltonian ([T]) and ignore terms containing 
nonzero momentum modes. We then have 



-J(<oVl/2^0 + h.C.) + ^^*l^ 



(3) 



where 'i'j^o annihilates a particle in wire i with zero mo- 
mentum {q = 0). The above Hamiltonian ^ can be 
mapped into a spin system via Schwinger representation: 
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Then the Hamiltonian Hq can be written as 

9 o2 



Ho 



'2iJ Sx 
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(4) 
(5) 
(6) 

(7) 



For a system of N atoms, Eqs. (H))-® map the bosons 
to the spin with the magnitude S = A^/2. We can easily 
find the Heisenbcrg equations of motion from the SU (2) 
algebra. Since we pick K large enough to be in the semi- 
classical limit and the zero momentum mode is highly 
occupied initially, we can safely ignore the effects of quan- 
tum fluctuations and treat Hq classically. Furthermore, 
since depiction is small, we can take N to be the total 



number of atoms in the system. Using energy and total 
spin conservation (particle number conservation in origi- 
nal fields), we can write down the equation of motion for 



2gEo 
L 



(8) 



where Eq is the total (conserved) energy of the system. 
Initially all atoms are in the first wire, so we have Sx = 
Sy = and Sz = N/2 as our initial conditions, which 
uniquely determine Eq. It is convenient to work with 
the dimensionless time r = gpot = gNt/L and rescale 
the spin Sz ^ n = Sz/S = {Ni - N2)/N so that it 
represents the scaled population difference between the 
two wires. Then the equation of motion reads (see also 
Ref. [H) 



din = - 4A^ - 
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where A = J/gpo = '^/m- The parameter A determines 
whether or not we have MQST. Eq. ^ is equivalent to 
the equation of motion for a single particle in the poten- 
tial 



V{n) 



4A^ 



n 

Y 



(10) 



and total energy E = [druf /2 + V{n). For A > l/VS 
the potential V{n) has a single minimum at n = (equal 
population in both wires), while for A < the min- 

imum at ?i = becomes a local maximum and two new 
minima appear at 



± n*, where n^, = ^/l — 8A^. 



(11) 



Nonzero n signifies population imbalance, so a necessary 
condition for MQST is A < But this condition is 

not sufficient for MQST, since a system with high enough 
energy E > V{0) = will overcome the potential barrier 
and not self-trap. We also need E < 0. For our initial 
conditions of complete imbalance (ri(0) — 1, i9tn(0) = 0) 
we have E = 2A^ — 1/8, so in order to have MQST we 
need to have 



A < 1/4. 



(12) 



In Figure [T] we show the effective single particle poten- 
tial, along with the energies, for A = 1 where there is no 
self trapping, and A = 1/5 where we have MQST. We 
can write down the exact solution for Eq. ^ in terms of 
a Jacobi elliptic function 



n{T) = cn ( 2At, ^ 



(13) 



Since we are interested in the self trapped regime, where 
A is small, we can expand Eq. (|13p to second order in A 
and ignore higher order terms 



n(r) ~ 1 - 8A 



2sin2 



~2~' 



(14) 
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FIG. 1: (Color online) Plot of V{n) for A = 1 (top, not self 
trapped) and A = 1/5 (bottom, self trapped), along with the 
energies E{X) for both cases. Note that we have the constraint 
— 1 < 71 < 1 from total number conservation. 

where n* is given by Eq. This approximation is very 
good for A as high as 1/5 (recall that we need A < 1/4 
in order to have MQST). In MQST we only have a par- 
tial tunneling of atoms between the two wires. Atoms 
oscillate back and forth between the wires, with a fre- 
quency n^, ~ 1 in dimensionless units (which corresponds 
to the frequency /i in original units) and the amplitude 
4A^. So in MQST we have high frequency and small am- 
plitude transfer of atoms between both wires. Whenever 
the effect of quantum fluctuations is not important, this 
situation persists and MQST is stable. As we described 
above these condensate oscillations couple to the non- 
condensate modes. Our next step will be to analyze the 
effect of this coupling on dynamics. 



B. Effective Hamiltonian for nonzero momentum 
modes 

We will now focus on finding the effective Hamiltonian 
for the nonzero momentum modes. We will employ the 



following strategy, as explained in the introduction of this 
section. We go back to the microscopic Hamiltonian and 
treat the zero momentum mode fields as classical time de- 
pendent sources coupled to nonzero momentum modes, 
their time dependent behavior given by the solution to 
the equations for zero momentum mode fields described 
above. Wc then keep contributions from nonzero momen- 
tum terms up to quadratic order, and solve the equations 
of motion for all nonzero momentum modes. The dynam- 
ics described by this effective Hamiltonian then should be 
valid as long as the occupation in the nonzero momen- 
tum modes remains small. Through this effective Hamil- 
tonian, we will find the mechanism of MQST decay and 
illustrate the importance of quantum fluctuations. We 
will then need to use the explicit solution for the zero 
momentum mode fields ^i,o(t) and '^2.o(t), which we 
can easily find from solving the classical system given by 
Eq. ([7]). For A < 1/5 to a very good approximation (up 
to terms of the order of A''') the solution is 




where (j)i^o is the initial phase in the zero momentum 
mode for wire i and N is the total number of atoms. 
Since the our initial state is the ground state of decoupled 
wires, these phases are random and in order to compute 
the expectation value of an arbitrary operator one needs 
to take average over them (technically this follows from 
the fact that in the classical limit the Wigner function 
representing initial Fock state reduces to the fixed am- 
plitude of 4* 1,2 and the random phase (boI. [31|). But in 
our particular case these phases are not important be- 
cause initially the second wire is empty so one can safely 
set them to zero. Note Eqs. (fTSj) and (fT6|) have a running 
phase which is typical of states in MQST. Also in the 
range of A < 1/5 we are interested in we can set = 1 
to a very good approximation. 

Next we will expand the Hamiltonian to the quadratic 
order in non-zero momentum modes. Doing this we will 
also ignore contributions to the Hamiltonian of order A^ 
and higher. This approximation is justified for small A, 
and small enough system size such that £\ = JL / n5,1^ 1 , 
where ^ is the coherence length in the full wire. The con- 
straint on system size is due to the fact that we need 
characteristic kinetic energies in the system (ex l/t"^) to 
be larger than the terms of order A^ ignored in the Hamil- 
tonian. This requirement on the system size not only 
simplifies the analysis, but (as we show later) is also cru- 
cial in order to obtain a sharp momentum distribution at 
intermediate times (r ~ 1/A). 

The effective Hamiltonian for the nonzero momentum 
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modes can be written as H = "^^^q Hq. As before, we 
will employ dimensionless units (or equivalently, all ener- 
gies in terms of /x, the natural inverse time unit). There 
are three contributions to Hq , which we will consider and 
discuss separately: 



m + m 



(17) 



where is the contribution from the full wire, is 
from the empty wire, and i/^ is the contribution from 
tunneling between the wires. For the full wire, we have 



for the coupling between the two wires, which is simply 

Hq' = -A(vI/t^vl/2^ + *I_^vI;2_^ + h.c.). (20) 

The explicit time dependence of Hq can be eliminated 
by a simple unitary transformation: 



Hq^UqHqUl+l{drUq)Ul 



where 



(21) 



Uq = exp lTY,{^>l^^.,,q + ^l_q^^^-q) , (22) 



^ Iq^ 1, — g 



■2ir 



-(2 + r,)(vi/t *i, + vi/t *,^_^), (18) 



where Tq = /2miJL, the scaled kinetic energy for mo- 
mentum q. Note that we have an explicit time de- 
pendence in this contribution coming from the running 
phase. The phase monotonically increases at a rate 1 in 
dimensionless time (with rate /i in original units). This 
is also the frequency with which the atoms tunnel be- 
tween wires in the self trapped regime. Physically, this 
contribution couples an external harmonic source with 
a frequency /i to creation and annihilation operators for 
particle pairs with momenta ±g. This acts like an ex- 
ternal energy pump, and is the cause of the dynamical 
instability that develops in the system due to resonance of 
the external frequency ^ with the natural system modes. 
For the empty wire, we have 



= Tq{-^\q^2q 



(19) 



Note that no explicit time dependence shows up in this 
contribution if we ignore terms of order A^, and that 
Eq. is just the Hamiltonian for free Bosons of mo- 
mentum q. There is a subtle issue in ID, where in the 
ground state of a dilute interacting Bosonic gas the in- 
teraction energy is always higher than the kinetic energy, 
since the interaction energy is proportional to density 
while kinetic energy is proportional to density squared. 
At short times the density in the empty wire is very small, 
so naively it seems that we wrongfully ignored the inter- 
action. However, this consideration is somewhat mis- 
leading for two reasons. First, we have included effects 
of interactions in the empty wire already: recall that 
we treated the zero momentum modes separately, and 
in this treatment we have included the effect of interac- 
tion in the zero momentum modes. As long as depletion 
from the quasicondensate is small, which it is for cases 
we consider, the interaction terms involving nonzero mo- 
mentum modes are negligible compared to their kinetic 
energy. Second, we are dealing with far from equilibrium 
systems and the momentum modes which are populated 
have very high energy of the order of the chemical po- 
tential in the full wire. While the particle density in the 
initially empty wire remains small the interaction energy 
is negligible compared to the kinetic energy and can be 
safely neglected. Finally we are left with the contribution 



where the second term in ((2T|) reflects the fact that Uq is 
time dependent. Physically this comes about because we 
are going to an accelerating (rotating) frame. The new 
Hamiltonian Hq no longer has any explicit time depen- 
dence. 

We then have for each of the contributions in turn 



-f(i + r,)(vi/t^*i, + *l 

H^ = {Tq-l){-^l^^2q 
Hq = -X{^\q^2q + ^Vq^2- 



(23) 
(24) 
(25) 



Note that the unitary (gauge) transormation simply re- 
sults in Tq — > — 1 i.e. subtracting the chemical po- 
tential of the full wire from the kinetic energy. The dif- 
ference in sign in the term T^ ± 1 between the full and 
empty wire has a simple physical interpretation: in order 
for a particle to tunnel from the empty wire to the full 
wire it must overcome the difference in interaction ener- 
gies in the wire. In other words, the chemical potential 
difference between the two wires acts as a bias (uniform) 
potential, with the full wire at a higher potential due to 
its interaction energy. 

Since Hq is quadratic, we can explicitly diagonalize it. 
Instead of bosonic fields ^'i „ and „ it is convenient to 

i.q 

work with the fields which represent density and phase 
fluctuations: 



_ ni,q + i(j)i^q 



V2 



(26) 



Since this is a canonical transformation the new fields sat- 
isfy standard coordinate-momentum like commutation 
relations: 



K,g, (l>j,q'] = iSq-q'S^j 



'^i,q ~ ^i,-qi 



4. 



(27) 



Their Fourier transforms ni[x) ~ l/T n^.^ exp(igx) 
and (f>i{x) = '^^4ii^qCxp{iqx) satisfy similar commuta- 
tion relations in real space 

n^{x)'^ = ni{x) 4>i{xy = 4>i{x). (28) 
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which are the correct relation for physical number and 
phase operators. 

To see that the fields ni^q and 0i,g indeed correspond 
to the correct number and phase operators for small de- 
pletion we observe that the Fourier component of the 
density 



(29) 



basis ■tpq 



'19 



V2 



(34) 
(35) 



where ipq and V'J obey the usual commutation relations 

[V'gj'/'g'] = ^qq'- terms of quasiparticles is given 
by 



is exactly proportional to ni_q. Since (j)i_q represents the 
conjugate field to then up to a constant factor 

{I/V2N) it is just the Fourier component of the physical 
phase (this can be checked by explicit calculation). The 
approximate equality above is only valid when the de- 
pletion is small and the terms with q' = 0, —q dominate 
the sum. Note that this argument cannot be applied to 
n2q and (j)2q, the fields in the second wire, since we have 
a small occupation in the zero momentum mode. n2q 
and (j)2q therefore do not correspond to physical number 
and phase fluctuations in the empty wire, but for sim- 
plicity we will use the same terminology in referring to 
them, since their properties are otherwise identical to the 
number and phase operators in the full wire. 

We can then write our Hamiltonian as Hq = + 
+ H-^ where 



H. 



^Tq{2 + Tq) (^^% + ^l^ij-q 



(36) 



Tq(piq(l)i^^q + {Tq + 2)ni,ni,_, (30) 

{Tq - 1) {n2qn2-q + (j)2q<l)2-q) (31) 
-A {niqn2.-q + n2qni^^q 

+4>lq4>2,-q + 4>2q4>l,-q) , (32) 



We can easily find the expectation values of the occupa- 
tion of momentum modes by real particles in the full and 
empty wires using the following relation 




The energy spectrum for the full wire ujf{(i) = 
yjTq{2 -\- Tq) \s just thc usual Bogoliubov spectrum for 
Bose gases close to thc ground state. In the original en- 
ergy units this spectrum reads 



where Vs — \f[ijrn is thc speed of sound in thc full wire. 
For small momenta topiq) — Vs\q\ and we have thc usual 
sound modes for a Bose gas close to zero temperature. 
Note however that we start in a high energy state with 
respect to the coupled Hamiltonian (interaction energy ~ 
fi) , and we need to take into account the correct spectrum 
for high momentum modes q ^ mvs as given by Eq. (). 

In order to treat perturbatively, we will write it 
in the tp, ^'2 basis, which diagonalizes and re- 
spectively. In this basis, iJ^ has terms with all possible 

quadratic combinations of tp,'^2,4'^ ,"^2- Consider that 
our initial state has very small depletion from the zero 
momentum modes, and therefore all contributions of 
with at least one annihilation operator have vanishingly 
small matrix elements compared with contributions of 
with two creation operators. So for our pcrturbative 
calculation to a very good approximation 



C. Dynamics of nonzero momentum modes 



iTqi2 + Tq)f' 



ii^-q'^lq+4'H^-q), (38) 



The Hamiltonian Hq in Eq. ((321) is quadratic and there- 
fore it is easy to analyze the dynamics in our system. 
Before we go on to solve the system Hq exactly let us 
gain some physical insight by doing a pcrturbative cal- 
culation for the system treating H^ as the perturba- 
tion. In the empty wire the Hamiltonian H^ is di- 
agonalized in the original particle basis (sec Eq. (|24p ). 
Its spectrum is a;_E(g) ~ Tq — 1 (in the original units 
^E{q) = 9^/(2m) — Apart from subtraction of the 
chemical potential due to our gauge transformation, this 
is just the energy spectrum for free bosons. 

Next, let us diagonalize the Hamiltonian for the full 
wire H^ ((50)) . This can be done in the quasiparticle 



It is clear from Eq. (|38p that the physical process that 
drives the depletion from thc zero momentum modes and 
destroys MQST is the production of quasiparticle pairs 
with opposite momenta ±g in each wire. At long times 
the dominant contribution to particle creation will come 
from the momentum mode satisfying the resonant condi- 
tion 



UE{qc) + C^F(gc) = {Tq^_ - 1) + ^Tq^{2+TqJ = (39) 



which has the solution 



2mfi 



(40) 
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In the original units Eq. ([39)) reads 



Q2_ 
2m 



with the solution 




{2mVg 



mVs 

71" 



(41) 



(42) 



This wavevector describes the characteristic momentum 
where the zero momentum mode particles scatter with 
the highest probability, possibly causing a breakdown of 
MQST. As we will see later Qc represents the momentum 
of the maximally unstable mode. 

Note that Eq. (|^T|) can be rewritten in the following 
way, giving a clear physical interpretation: 



q2 / q2 
— +Vs\q\\ 1 + — 



fj.. 



(43) 



The right hand side of Eq. ((43|) is just the sum of the 
energies of two quasiparticles of momentum with one 
quasiparticle in the full wire and the other in the empty 
wire, while the left hand side is just the frequency of 
oscillations of the zero momentum mode, which serves 
the external source term. Through the resonant coupling 
particles are scattered out of the zero-momentum mode 
to the momenta close to Qc- Because of the bosonic na- 
ture of quasi-particle, we find that the population of the 
mode with q = qc increases exponentially with the rate 
2 J/>/3. If the number of unstable modes is small we can 
expect a very sharp momentum distribution forming in 
both empty and full wires at time scales of the order 1/ J. 

Now we outline the exact solution for the dynamics 
of the system in the linearized regime. We can easily 
find the Heiscnberg equations of motion for the nonzero 
momentum modes using Eqs. (|30 p -(|32 |) 



drn2q 
dr<l>2q 



Tq(t)lq 
-{Tq- 
(Tq- 
-(Tq- 



2)niq + Xn2q 
)4>2q - Mlq 

1)^2, -I- Aril,. 



(44) 
(45) 
(46) 
(47) 



Assuming the solutions have a time dependence ^ 
Q-iuj(q)T ^ we find two pairs of eigenfrcquencics for each 
momentum: -iiu}±{q) (measured in units of /i): 



:(9)=r2+A 



1 1 

2±2 



(4r, - 1)= 



4A2(4r|-l). 

(48) 

Note that in the self trapped regime (0 < A < 1/4) the 
term under the radical in Eq. (|48p can be negative. Thus, 
there is always an interval in Tq where lo± is complex. 

In order to get a better physical understanding of the 
expression in Eq. (|48p let us first analyze the trivial limit 
A = where Eq. ([48)1 reduces to 



UJ 



(49) 




FIG. 2: (Color online) Plot of as a function of Tq at A = 0. 
Solid red line corresponds to uj'L and dashed green line cor- 
responds to Note that the two curves cross at Tq = 1/4, 
which is the resonance point (see Eq. (|40}) where the fre- 
quency of oscillations in the zero momentum modes matches 
the energy required to create particle pairs of characteristic 
momenta ±gc (see Eq. (|42|l ') in each wire . 



A,=1/8 




FIG. 3: (Color online) Same as in Fig.[lbut for A = 1/8. In 
this case uj± acquires an imaginary part in the region around 
the resonance point (see Eq. (|40|l ). so we have a dynamically 
unstable set of momenta in this region. The extent of the 
unstable region is ^/3A/2 while the maximum peak for the 
imaginary part of ijj± is A/\/3. The momentum modes falling 
within the unstable region then become populated at an ex- 
ponential rate, in a time scale ~ 1/A in dimensionless units 
( 1/ J in original units) . 



The two branches of the spectrum Lo\{q) are plotted in 
Figured] The two curves cross at = 1/4, the res- 
onance point we found in the perturbative calculation 
(see Eq. HO]) . The interpretation of uj± at A = here is 
simple. For < < 1/4, = {Tq - 1)^ corresponds 
to the spectrum of the free Bose gas in the empty wire 
shifted to the right by an amount equal to the chemi- 
cal potential difference between the wires. Likewise, for 
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< Tg < 1/4, a;_ is the spectrum of the interacting Bose 
gas in the full wire, and for Tq > 1/4, lo^ and uj- switch 
roles. For A = 0, the two wires do not interact with 
each other and their spectra are completely independent. 
For small nonzero A, however, the crossing point in Fig- 
ure [2] gives the resonance condition we found from the 
perturbativc calculation. It is easy to see that this is the 
channel for MQST decay. 

Next we take A small but nonzero. We plot uj± for 
A = 1/8 in Figure [3l In this case uj± acquires an 
imaginary part. This signals the development of a dy- 
namical instability in the region where U!± is complex. 
The dynamically unstable region is centered around the 
resonance point found in our perturbativc calculation 
{Tq = 1/4), and the extent of the region is ~ \/3X/2. In 
terms of momentum, the dynamically unstable region is 
centered around qc = mVs/V2 (where the speed of sound 
Vs = and the extent of the region is ^/SXqc- From the 
exact spectrum of Hq wc then find that there is an unsta- 
ble region of finite extent ^ A around the resonance point, 
instead of the single characteristic mode found from the 
perturbativc calculation. It then follows that the popula- 
tion of all momentum modes that fall within the unstable 
region, will exponentially grow in time. The time scale in 
which this happens is given by the inverse of the typical 
value of the imaginary part of uj± , which is of the order 
1/A (1/J in original energy units). The maximum rate 
of increase possible is determined from the peak of the 
imaginary part of lj±, which is 2J/a/3. This is precisely 
what we obtained from the perturbativc Fermi Golden 
Rule calculation. 



D. Second dynamically unstable region 

Let us briefly discuss other regions of interest in the 
system. From Eq. we find that there is an addi- 
tional instability occurring around 2qc, which is twice 
the resonant momentum. However, this time only uj-, 
and not becomes complex. We then have a second 
region of dynamical instability in the system. It is easy 
to see that there are no additional unstable regions in 
the system. This instability emerges because cj^(2gc) be- 
comes negative for any small positive A. This instability, 
however, is much weaker since the maximum imaginary 
part of the frequency is A^/3. Since the main instability 
around qc develops on times of the order of 1/A the ef- 
fect of the second instability on dynamics is usually very 
small. Nevertheless, as we will see from full nonlinear 
numerical simulations, this instability is not just an ar- 
tifact of the effective model and the corresponding mode 
is indeed excited in the system. In Fig. [?] we plot imagi- 
nary part of uj-iq) for A close to 1/4, i.e. to the onset of 
self-trapping. It is clear this instability also has a smaller 
width (by a factor A/3\/3) compared to the first insta- 
bility. For smaller values of A the effect of the second 
instability is even weaker. 

Let us briefly mention that the physical process asso- 



ciated with this instability can be understood by looking 
into perturbation theory. Since the instability here is 
on the order of A smaller than the one associated with 
the first unstable region, we expect the effect to show up 
in second order perturbation theory where we have con- 
tributions of order A^. This is indeed the case, and we 
find there is a resonance here similar to what occurs in 
first order. In second order, we have a contribution pro- 
portional to the matrix elements of ^'J g'!/'g'0j^2 -g ^^"^ 
likewise for q —q. This contribution has the net ef- 
fect of producing particle pairs of opposite momenta in 
the empty wire, where the resonance condition for these 
terms is now given by 2{Tq — 1) = instead of Eq. ((M)) . 
The characteristic mode then satisfies 2{q^/2m — /i) = 
which gives ns q ~ 2qc. This process can be interpreted 
as two atoms in the full wire zero momentum mode (with 
interaction energy fj, per atom) tunneling into the empty 
wire, where all of their interaction energy is converted 
into kinetic energy. There is a similar process in the full 
wire as well, where two atoms in the zero momentum 
mode are transfered to the 2qc mode and in the pro- 
cess all their interaction energy is converted into kinetic 
energy. One can show that to all orders of perturbation 
theory only and 2qc are resonantly populated, which is 
consistent with the exact solution of the effective model. 



E. Effect of system size 

The finite system size plays a crucial role in our anal- 
ysis. First, as we noted earlier, we are relying on having 
a large condensate fraction in the initially full wire. This 
is only possible in one dimension if the system size is 
finite and small enough (see Eq. Additionally fi- 

nite system size leads to quantization of allowed momen- 
tum modes. When we analyzed the spectrum of unstable 
modes above we found that there are always unstable re- 
gions corresponding to imaginary frequencies. However, 
having these unstable regions does not imply that a finite 
size system actually has unstable modes and always de- 
cays from MQST. The dynamical instability rather pro- 
vides us with a possible route of MQST decay. The more 
accurate statement is that weakly interacting coupled ID 
systems that initially exhibit MQST will undergo MQST 
decay through the resonant mechanism described above 
as long as there are actual momentum modes in the sys- 
tem that fall within the dynamically unstable region. 
The set of allowed modes is controlled by the system 
size, so the system size will determine whether or not we 
have decay from MQST. 

In the following it is easier to work with the system 
size in terms of the coherence length [l = L/^). It is 
also convenient to measure the momentum in units of 
1/^, so that the allowed momentum modes are given by 
qn = 27rn/£ with an integer n (33 |. The characteristic 
momentum mode that is resonant with the zero momen- 
tum mode oscillations is given by q^. = 1/^/2 in these 
units. The condition for MQST decay is simply stated 
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A = 1/4 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

q (units of 1/^) 



FIG. 4: (Color online) Plot of (in units of vs q (in units of 1/^). Unlike tj+ where there is only one region where it 

acquires an imaginary part, for uj- we have two such regions. The first region is exactly the same as in u}+, centered around 
q = qc = l/\/2^. The second unstable region occurs only for uj- around q = 2qc. Both the extent of the region of instability and 
the maximum of the imaginary part of oj in the second region are of order A times smaller than the corresponding quantities 
in the first region, meaning that the instability in the second region is suppressed compared to the one in the first region if A 
is small. 



that for a given system size i, there exists at least one 
value of the integer n such that 



27m 



(50) 



If this condition is not satisfied for any n then MQST is 
stable even for a ID system. 

The requirement Eq. [SD] implies that there is a mini- 
mum system size ^min below which MQST is always sta- 
ble: 



27r 



qc{l 



2 J 



(51) 



For small A to a very good approximation wc find that 
this is equivalent to the minimum system size (in original 
units) imin = 2n\/2£_. For all systems smaller than ^min 
MQST is always stable and the effect of quantum fluctua- 
tions can be ignored. We then have a situation identical 
to the Joscphson junction where MQST was originally 
studied, where only the dynamics in a single momentum 
mode is important in each wire. 

We also note that there exists a system size ^thresh 
above which MQST is always unstable and the dynam- 
ical instability is always present. This corresponds to 



the situation whenc there is always at least one solu- 
tion satisying Eq. ([SO)) . The condition that determines 
the threshold length is that the spacing between adjacent 
momentum modes (that is 2tt/£) is smaller than the size 
of the dynamically unstable region in momentum space 
(see Figure S]). The threshold length is given by 



^thrc 



2tt 



(52) 



or in original units it is just itiuosh = T'y f^' 
threshold length is of the order of ^/A. The smaller 
we make the tunneling between the wires, the larger the 
threshold length will be. 

Finally, there is a range of intermediate system sizes 
between the minimum size ^min below which MQST is 
always stable, and the threshold size ^thresh above which 
MQST is always unstable. For system sizes £niin < i < 
^thresh, WC find that by changing the system size we al- 
ternate between intervals where MQST is stable and in- 
tervals where MQST is unstable. Moreover in unstable 
intervals there is always a single unstable mode. It is 
easy to see that unstable intervals in £ (for intermediate 
system sizes) are centered around 



27rn 



(53) 
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where n is an integer. If we define nmax as the largest 
value of n for which ^„ < i'thresh is satisfied, then we need 
to have rimax < "^/f^- We also note that the size of the 
interval in t. where MQST is unstable increases linearly 
with n so that the n-th unstable region corresponds to 
the following interval for t. 

i^<£<i^. (54) 

For n > riinax, the size of the unstable regions become so 
large that they overlap with each other and the system 
is always unstable. 

Because in the intermediate regime i^-m < i < ^thresh 
there is at most one unstable mode, we should be able 
to observe a very sharp momentum distribution develop 
in unstable regions. Furthermore we can control how 
quickly the system decays from MQST; specifically we 
can set the system size such that the rate of decay has 
a maximum. This happens if the system parameters 
are such that the unstable mode perfectly matches the 
resonance condition. Clearly by decreasing the tunnel- 
ing A we can make resonance conditions sharper and 
the intermediate region larger. To achieve the maximal 
growth rate in the population of the resonant mode (and 
hence the maximum number of particles transfered to 
this mode) we need to choose £ = = ^nn/qc. Note 
that by picking such in we are also guaranteed to get a 
momentum state that falls within the second unstable re- 
gion q ~ 2qc in momentum space (see Figure Recall, 
however, that the rate of growth here is of the order of A 
smaller than the corresponding one in the first unstable 
region, so we may neglect its effect if A is small enough, 
which it is in cases where we initially have MQST. Let 
us finally mention that if the system is in the intermedi- 
ate size regime, then according to Eq. ((M]) (see Figure 0] 
as well) one can easily tune the system from stable to 
unstable MQST by tuning A (or equivalently J). 

The analysis we presented so far neglects the interac- 
tions between different modes. In particular, it neglects 
the damping of the zero momentum mode oscillations, 
which is inevitable due to creating high energy excita- 
tions in the system. This damping will destroy resonant 
population of the unstable modes and thus will lead to 
eventual saturation of the exponential growth. At longer 
times one can expect that the energy stored in the unsta- 
ble momentum mode will be redistributed among other 
modes and the system will eventually thermalize (or at 
least reach some steady state). The two important ques- 
tions are (i) What is the maximum achievable occupation 
of the resonant modes? and (ii) What is the time scale 
at which the sharp distribution disappears and the sys- 
tem thermalizes? Later when we do numerics, we will see 
that higher n decreases the thermalization time for the 
system and decreases the maximum occupancy of the res- 
onant modes. Therefore by picking smaller system sizes 
we can delay thermalization. If we go to the extreme 
case and pick the smallest system size, so that only the 



first excited momentum mode falls within the dynami- 
cally unstable region, we do not observe thermalization 
in our numerics for the time scales considered. In particu- 
lar, in this regime we do not see equilibration of particles 
between the two wires, and the system seems to remain 
self trapped for very long times. In this extreme case 
the sharp momentum distribution is stable, and the sys- 
tem seems to reach a steady state in population of the 
nonzero momentum modes. We do observe, however, a 
slow leak of atoms from the zero momentum mode of the 
full wire to the zero momentum mode of the empty wire, 
as seen in Figs. [S] and [H] Only the zero and first excited 
momentum modes are important in the dynamics. This 
is in contrast to larger system sizes (with modes in the 
unstable region), where the sharp momentum distribu- 
tion disappears, the population difference between wires 
disappears, and the system thermalizes at much shorter 
times scales. 

In general thermalization times decrease as the sys- 
tem size increases. If we pick system sizes above ^thresh, 
MQST always decays, and if £ is large enough we will 
always have multiple modes fall within the dynamically 
unstable region. If enough modes fall within this region 
we will no longer have a sharp momentum distribution 
in intermediate times, and in fact will not even be able 
to discern that there is MQST in the system. The time 
scale for MQST decay is on the order of J if only one, 
or at most a few modes fall within the unstable region. 
This time scale decreases as we make £ larger until even- 
tually it becomes on the order of /i, the frequency of zero 
momentum mode oscillations. At this point, we can no 
longer claim that we had MQST in the first place. So in 
the thermodynamic limit there is never any MQST in ID 
systems. 

Let us note that while our analysis directly applies to 
uniform systems, one can expect that qualitatively the 
results will be the same in a parabolic trap. In this case 
the noncondensate eigenmodes will not be plane waves 
but wavepackets oscillating with characteristic momen- 
tum, which will be quantized. If the resonant conditions 
arc satisfied we expect MQST decay while if they are not, 
MQST wiU be stable. 



III. NUMERICS 

In this section we will present the numerical results. 
Our system is described by the microscopic Hamiltonian 
P]). We will work with weakly interacting systems. In 
this regime the initial condensate fraction is large and the 
resonant mode population is most pronounced. Also for 
weak interaction we can reliably use truncated Wigner 
approximation (TWA) method to simulate the dynam- 
ics (see Refs. @, IH, [H, [H H, for details of this 
method). TWA gives the leading order in expansion of 
quantum dynamics in quantum fluctuations around the 
classical (Gross-Pitaevskii) limit. The classical limit can 
be strongly interacting in the sense that classical dynam- 
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ics is intrinsically nonlinear. It can be shown that in 
TWA the classical equations of motion do not change 
but the initial conditions for classical fields become fluc- 
tuating and distributed according to the Wigner func- 
tion. These fluctuations are crucial for our problem since 
without them non-zero momentum modes are never ex- 
cited. We emphasize that we use the term "classical" for 
Gross-Pitaevskii equations rather than mean-field com- 
monly used in literature. It is important to realize that 
in TWA each classical trajectory does not describe any 
mean-field, it rather describes characteristics along which 
the Wigner function is approximately conserved. Only 
in the classical limit where there are no fluctuations 
and all particles are described by a single condensate 
mode Gross-Pitaevskii equations describe mean-field dy- 
namics. Effectively the Wigner function in TWA gives 
the quasiprobability distribution for classical fields deter- 
mined by the initial state of the system, pure or mixed. 
At higher orders, there will be also corrections to the 
classical trajectories, but in the regime of large K we are 
interested in, those are unimportant p?!. [29|. 

The implementation of TWA is straightforward. The 
dynamics of the system is described by Gross-Pitaevskii 
equations of motion for the following Hamiltonian 

^Jo V2m 2 

^ J / da; {■^\-^2 + *2*i) , (55) 
Jo 

where we simply replaced the quantum operators in 
Eq. Il]) with classical fields. We note that in principle 
one has to use the Weyl symbol of the Hamiltonian in 
Eq. (j55p . However, for spatially uniform quartic inter- 
actions the difference between ([55]) and the Weyl sym- 
bol is given by a constant term proportional to the to- 
tal number of particles. In equations of motion it yields 
an overall phase, which is not important. The Gross- 
Pitaevskii equations corresponding to this Hamiltonian 
are obtained in the standard way: 



(56) 



Note that the same equation can be obtained by intro- 
ducing an analogue of Poisson brackets for the coherent 
states and writing standard classical Hamilton equations 
of motion [2^ . 

To implement TWA we need to solve equations ((56|) 
supplemented by random initial conditions distributed 
according to the Wigner transform of the initial density 
matrix. Since we are dealing with an interacting system 
in the initial state, finding the Wigner transform is not 
trivial by itself. One possibility to overcome this diffi- 
culty is to find the initial state within the Bogoliubov's 
approximation and then find the Wigner transform of 
the corresponding wave function psj . We will employ a 
different route. Namely we start with a noninteracting 
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FIG. 5: (Color online) Momentum distribution in the empty 
wire for L — 2n/qc, the smallest system size for which MQST 
decay can be observed, and A = 0.075. The two sets of 
curves correspond to the Luttinger parameters K/n — 50 and 
K/ir = 100 respectively. Short time small amplitude oscilla- 
tions of g = modes correspond to MQST. At intermediate 
times resonant modes ai q = Qc are populated and the mo- 
mentum distribution is sharply peaked at q — ±gc reaching 
the maximum population of about 8% of particles per mode. 
The maximum is bigger for weaker interactions (larger K). 
After the saturation of the resonant mode we observe very 
slow leak of atoms into the empty wire in the zero momen- 
tum mode with a large population imbalance persisting for 
very long times. Note that only few modes {q = 0, ±gc, ±2gc) 
are involved in the dynamics. 



full wire (g — 0) where it is straightforward to find the 
Wigner transform of the ground state. Then we slowly 
(adiabatically) turn on interactions in the system (keep- 
ing the wires decoupled) until we reach full interaction 
strength. Since quantum fluctuations remain weak in 
the whole range of g this procedure gives us the correct 
initial state for the interacting system (we checked that 
this procedure is very accurate on a small lattice system 
comparing such adiabatic ramp with exact simulations). 
We note that in this setup it is actually very hard to 
achieve adiabatic limit in large ID systems |32l |. How- 
ever, this does not represent a real issue for us since we 
confine ourselves to mesoscopic system sizes. Once the 
desired interaction strength is reached we suddenly turn 
on coupling between the wires and analyse the subse- 
quent dynamics. 

We consider weakly interacting systems with K/tt = 
50, 100 so that we are certainly in the semiclassical limit. 
Since K/tt = po^ we see that this corresponds to the 
regime of roughly 100 particles within the healing length. 
The minimal length of the wire required to see decay of 
MQST is Lmin = 27rV2^ « 9^ (see Eq. i^) thus we see 
that this parameter regime requires of the order of few 
hundred particles in the full wire. The effect will persist 
for smaller values of K/tt ^ 10 so it can be observed 
with smaller number of particles per wire, though for 
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Full wire momentum distribution (A,=0.075) 
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FIG. 6: (Color online) Same as Figure [5] but for the full wire. 
Again, only the lowest modes are involved in the dynamics 
{q = 0, ±gc, ±2gc) and we find no equilibration in these time 
scales. 



smaller K the momentum distribution becomes not so 
sharp. We will choose small couplings A = 0.075, 0.05 so 
that the two relevant time scales in the system (1/J and 
l//i) are well separated and we can clearly distinguish the 
MQST part of the dynamics at short times from the sharp 
momentum distribution that develops at times ^ 1/J. 
Finally, we will work mostly with system sizes i = L/S, ~ 
^1, iio < 4hresh (sec Eqs. ([52]) ■ (|53|) ) ■ so that we are in the 
regime with exactly one momentum mode falling within 
the dynamically unstable region. We also choose system 
sizes that maximize the imaginary part of the frequency 
of the unstable mode. This gives us the largest transfer 
of particles to these modes. 

Let us first discuss in more detail what happens for 
L = Li ~ 2n/qc. We have plotted the results for mo- 
mentum mode occupation for all momentum modes for 
times up to i = 225/^ in Figure [5] for the empty wire, 
and Figure [6] for the full wire. We choose A = 0.075 and 
two different Luttinger liquid parameters K/n ~ 50 and 
100. At short times wc clearly observe MQST, which 
is characterized by the high frequency, small amplitude 
oscillations in the zero momentum mode. The behav- 
ior at short times is similar for both K/tt = 50 and 
K/tt = 100, but note that the damping of oscillations is 
stronger for smaller K, which is not surprising since this 
corresponds to stronger quantum fluctuations. Then at 
times ~ 1/J, we find a large transfer of particles to the 
characteristic momentum modes qc, which for the given 
system size is just the first excited momentum mode in 
the system. A larger fraction of particles is transferred for 
larger K , which is expected since larger interactions pre- 
clude a large population from developing at any mode. 
We also clearly observe transfer of particles to the 2qc 
mode. Recall from our previous discussion that there is 
also an unstable region around 2qc, but the rate of trans- 
fer there is of order 1/A smaller. All other momentum 



modes remain virtually unpopulated at all times shown, 
so the dynamics is confined only to these first few modes, 
and we find that this greatly delays the onset of thermal- 
ization (see Figure [H]). Note that at later times the occu- 
pation in qc saturates at about 8% of the total number 
of particles, while the population in the zero momentum 
mode monotonically increases and eventually surpasses 
the population in the characteristic modes. The system 
is thermalizing, albeit very slowly. 

Let us now focus on a larger system size L = Liq = 
207r/gc. This system size is such that the tenth excited 
momentum mode falls within the unstable region. We 
pick K/tt — 100 and A = 0.05. This choice of parameters 
allows us to delay thermalization so that we can clearly 
distinguish what happens in different time scales in the 
system. The results for the empty wire momentum oc- 
cupation are shown in Figure [71 (a 3D plot is shown in 
Figure lT^ . At short times we once again observe the zero 
momentum mode oscillations present in MQST. We ob- 
serve the transfer of particles to the characteristic mode 
qc at intermediate time scales ^ 1/J. Note that when 
the population at qc peaks, the momentum distribution 
is very sharply peaked around ±qc as we see in Figure [S] 
The system thermalizes at later times. In Figure [5] we 
compare the momentum distribution of the empty wire 
at different times. At t = 7r//i, which corresponds to a 
half period in MQST oscillations, we find a small number 
of particles at zero momentum while all other modes are 
virtually unoccupied. At intermediate times (t — 160//i, 
which is of the order of 1/J) we find the momentum dis- 
tribution very sharply peaked around ±qc (with a small 
peak at g = 0), while all other modes remain unoccupied. 
At late times {t = 400/ju) the system shows a thermal 
momentum distribution centered around q = 0. Note 
that the thermalization time scale is much smaller here 
than when L = Li. Also the height of the peak when 
the population at qc is maximum is higher than the max- 
imum height of the thermalized momentum distribution. 
This fact taken with the fact the the distribution is very 
sharp when the population at qc peaks should make the 
intermediate sharp momentum distribution state easily 
distinguishable from the final thermalized state. For the 
latest times shown t ~ AOO/fi, wc compare the momen- 
tum distribution of both wires in Figure [HI Note that 
they arc virtually indistinguishable at this point so the 
system fully thermalized. 

There is an important difference between dynamics in 
short systems with £ ~ £i and longer systems with £ = 
where 2 < n < rithresh (compare Figs. [S] and [T]). While 
in both cases there is only one unstable mode which gets 
populated at intermediate times ^ 1/J the longer time 
dynamics is clearly different. In both cases the popula- 
tion of the resonant mode q ~ qc saturates at a value of 
the order of a few percent of the total population. At 
later times in the shorter wire the population remains 
confined essentially to the condensate (zero momentum) 
and q = ±qc momentum modes, while in longer wires 
other momentum modes get excited and the system thcr- 
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Empty wire momentum distribution (L=L^q) 
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FIG. 7: (Color online) Momentum distribution in the empty 
wire where we have L = 20n/qc, X — 0.05, and K/n — 100. At 
short times there is clear MQST then at times ~ 1/ J we find a 
very sharp momentum distribution as particles are transferred 
to the characteristic momentum modes. The system size is 
picked such that the characteristic mode corresponds to the 
tenth excited momentum mode in the system. At later times, 
particles are distributed in many momentum modes and the 
system thermalizes. This thermalization is accompanied by 
MQST decay. 
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FIG. 8: (Color online) Momentum distribution in the empty 
wire at several different times. We use the same parameters 
as in Figure [7] At time r = tt (corresponding to half the 
period of MQST dynamics) only the zero momentum mode is 
populated. At time r — 160 we find a very sharp momentum 
distribution centered around characteristic modes ±qc- This 
time corresponds to the peak population in dzq^ modes. Note 
that at this time the population in ±gc by far eclipses the pop- 
ulation of the zero momentum mode. Finally at time t — 400 
we observe a very broad thermal-like momentum distribution 
centered around q — 0. 



malizes. In order to observe sharp momentum distribu- 
tion in this case it is necessary that the thermalization 
time is longer than 1/ J, the time scale at which instabil- 
ity develops. We find that in order to ensure this is the 
case the following condition has to be satisfied: 



1 



Ln 



1/2 A' 



< 



4J2 

1^ 



(57) 



The right side of Eq. [57] is just the amplitude of oscilla- 
tions of the q = mode particle tunneling at short times 
(while MQST is present), and the left hand side gives us 
roughly the depletion of particles from the zero momen- 
tum mode in the full wire. The basic idea behind Eq. ((57|) 
is that the amplitude of MQST zero momentum mode 
oscillations should be easily distinguishable from the de- 
pletion background due to interactions in order to have 
well defined MQST dynamics at short times. We veri- 
fied numerically that if the inequality in Eq. (j57p is not 
satisfied, the momentum distribution is broader and the 
system thermalizes before the peaks can fully develop, 
making them harder to observe. Thus for A = 0.075, the 
RHS of Eq. ([FT)) is 0.025 so wc can have a few percent de- 
pletion at most in the full wire. If we increase A then we 
can tolerate larger depletion. However, A can not be too 
large otherwise oscillations of the condensate mode are 
unharmonic and the instability is suppressed. In practice 
we find that for A > 0.1 the system does not produce as 
sharp a momentum distribution. Even though the peaks 
are still distinguishable, their amplitude is smaller and 
they acquire some width compared to the case of smaller 
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FIG. 9; (Color online) Plot of the momentum distributions 
in both wires at time r = 400. Parameters are the same as 
in Figure [7] Note that the two curves are practically indistin- 
guishable. The system has thermalized at this point. 



A satisfying Eq. (|57|) . 

For completeness we show the momentum distribution 
for system sizes smaller than the minimum required for 
MQST decay in Figure [TUl (specifically L = tt/qc, which 
is half the minimum size required for MQST decay) and 
show that MQST is indeed stable in such a system. Note 
that even though the system remains self trapped, the 
zero momentum mode oscillations arc damped, which 
docs not happen in a single Josephson junction. The 
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FIG. 10: (Color online) Plot of zero momentum occupation 
for both full and empty wires for L = tt/qc. This is half the 
smallest system size necessary to have MQST decay. Note 
that the system remains self trapped, but the MQST oscilla- 
tions in the zero momentum modes are damped, unlike the 
usual case of MQST in a Bosonic Josephson junction. 
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FIG. 11: (Color online) Plot of zero momentum occupation 
for both full and empty wires for L — Stt/qc. According to our 
previous analysis this system should exhibit no MQST decay. 
This system size falls in an interval of MQST stability that is 
sandwiched between system size intervals where MQST is un- 
stable, which shows it is possible to tune the system size to al- 
ternate between regions of MQST stability and MQST decay. 
Note that the system remains self trapped, but the MQST 
oscillations in the zero momentum modes are damped, unlike 
the usual case of MQST in a Bosonic Josephson junction. 



source of damping likely originates from coupling to non- 
resonant modes, which are always present in ID (note 
that this system size still exceeds interparticle distance 
by two orders of magnitude). We also show that MQST 
is stable for L = Stt/qc in Figure [TT] This system size is 
exactly in the middle between the first and second insta- 
bilities. This shows that it is possible to tune the system 
size such that we alternate between MQST stability and 
MQST decay. Note that there is no observable thermal- 
ization in the system for very long times. Finally let 



us mention that for L ^ ^thresh (see Eq. ([52)) 1 where 
there are many unstable modes, the system thermalizes 
so quickly that there is no MQST dynamics present at 
all, and a sharp momentum distribution never develops. 

Consider now the very interesting possibility of decou- 
pling the wires immediately after obtaining the sharp mo- 
mentum distribution in the empty wire. Consider once 
again the parameters used in Figure [71 where we have 
L = 207r/(7c, K/tt = 100, and A = 0.05. Note that 
the population in the characteristic momentum modes 
peak around t ~ 140/^. To preserve this highly nonequi- 
librium momentum distribution at this moment we can 
suddenly decouple the wires. The resulting dynamics is 
shown in Figure 1121 Clearly the momentum distribu- 
tion remains peaked around the characteristic modes Qc, 
which is very different from the dynamics shown in Fig- 
ure [71 where the wires remain coupled. The system does 
not equilibrate in the usual sense. This is not surprising 
since after decoupling the wires we are left with a ID in- 
teracting Bose gas of the Lieb-Liniger type (see Ref [l^l ) , 
which is an integrable system solvable by Bethe ansatz. 
The integrability of the system precludes thermalization 
and as a consequence we observe that the sharp momen- 
tum distribution is robust. Through the procedure con- 
sidered here, we have found a way to construct a ID sys- 
tem with a stable sharp momentum distribution around 
nonzero characteristic momentum modes, which may be 
used itself in interesting applications. Note that in the 
empty wire almost all atoms occupy ±qc modes, which is 
a highly nonequilibrium state akin to having large coun- 
terpropagating currents in the same wire. In a usual 
nonintegrable system one expects that collisions between 
counterpropagating particles would thermalize the sys- 
tem and the large currents would dissipate, but integra- 
bility in the system considered here protects these cur- 
rents and we observe no current decay (see also Ref. [13 )■ 

Before we conclude this section, let us make a few com- 
ments on the regime of applicability of our results. We 
focused mainly in the case of fairly large K/tt ~ 50, 100 
since in these cases one has the very interesting situ- 
ation where we have a very large transfer of atoms to 
the characteristic modes and therefore a very sharp mo- 
mentum distribution develops before thermalization sets 
in, which makes this choice of parameters very interest- 
ing from a theory point of view. However, qualitatively 
our conclusions remain robust to much stronger interac- 
tions as long as Eq. ([2]) is satisfied. For smaller K we 
observe smaller transfer of particles to the characteris- 
tic modes and thus not so sharp momentum distribu- 
tion. We find that for K/tt > 10 the intermediate time 
momentum distribution is easily discernible from the fi- 
nal thermal distribution. For comparison, let us mention 
that for K/tt = 20 and £ — 207r/gc (so the characteris- 
tic mode is the tenth excited momentum mode), we get 
about 1.6% of atoms transferred to each characteristic 
mode (i.e. ~ 3% for both ±qc), compared to 4% for 
each characteristic mode for K/tt = 100, and the distri- 
bution remains fairly sharp for this case. For K/tt = 10, 
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FIG. 12: (Color online) Momentum distribution in the empty 
wire where we have L = 207r/gc, A — 0.05, and K/n — 100. 
These are the same parameters used in Figure [71 except now 
we decouple the wires att — 140/ fi, which is roughly when the 
population in the characteristic modes peaks. Note that the 
momentum distribution remains highly peaked around the 
characteristic momentum modes. Since the system is now 
integrable, we do not observe thermalization. 

however, the percentage transferred to each characteris- 
tic mode goes down to 0.6% and the distribution becomes 
noticeably broader. 

At the same time the nonmonotonic dependence of 
MQST decay on system size remains very robust as long 
as K/tt > 10. Thus for K/tt > 20, the system size de- 
pendence is practically indistinguishable from what hap- 
pens at larger K (i.e. K/tt = 100, see Figures [TU] and 
llip . In particular, we still see essentially no decay from 
MQST if the system size is such that there are no un- 
stable modes present, and if there are unstable modes 
present in the system we see decay of MQST in a similar 
fashion to what happens for larger K. For K/tt ~ 10, we 
still observe nonmonotonic dependence of MQST decay 
with system size, however unlike the situation with larger 
K (see Figs. [TUl [TT]) there is some small yet noticeable 
decay of MQST even when the resonance condition is not 
satisfied. For K/tt ~ 5, we find that the nonmonotonic 
behavior of MQST decay with system size is no longer 
present. 

IV. DISCUSSION 

In this paper we studied the breakdown of MQST for 
two coupled ID interacting Bose gases with a large pop- 
ulation imbalance. There are certain constraints that 
the system must meet in order for it to initially exhibit 
self trapping. In particular, the system size should be 
mesoscopic such that the initial depletion from the qua- 
sicondensate modes remains relatively small. Also the 
interactions in the system should be weak corresponding 
to large Luttinger liquid parameter K > 30. If these con- 
ditions are satisfied the initial self-trapping occurs if the 



ratio of the tunneling coupling J to the chemical poten- 
tial in the full wire /i is sufficiently small: A = J//J, < 1/4. 

Under these conditions the system exhibits MQST at 
short times. MQST is characterized by small amplitude 
(order A^) tunneling of zero momentum mode (conden- 
sate) particles between two ID wires. The frequency of 
these characteristic MQST zero momentum mode oscil- 
lations is /i, the chemical potential difference between 
wires. At intermediate times (of order 1/J, where J is 
the tunneling energy) we find that MQST can decay. The 
mechanism of this decay is the resonant coupling between 
the oscillating zero momentum mode fields and creation 
of quasiparticle pairs of characteristic momenta ±qc with 
the energy matching the oscillation frequency. These ex- 
cited modes have a very high momentum on the order 
of the inverse healing length in the full wire (or equiva- 
lently with the energy of the order of fi) . Through this 
mechanism, the system can develop a very sharp momen- 
tum distribution around ±gc at intermediate times. At 
later times, the system generally reaches a steady state 
at which point the momentum distribution in the two 
wires is identical. 

Since we are working with mesoscopic systems, we find 
that the system size puts a very strong constraint on 
whether or not we observe MQST decay, due to the fact 
that we have a discrete set of allowed momentum modes. 
Simply put, if there arc no modes that satisfy the res- 
onance condition (within a tolerance of order J), then 
MQST is stable even in ID systems and the population 
imbalance survives for very long times. A corollary to 
this is that there is a minimum system size Lmin — 27r/gc 
below which MQST is always stable, and a system size 
^thresh — 2tt / ^/SXqc abovc which one always has MQST 
decay. For intermediate system sizes between these two, 
in order to have MQST decay we must confine ourselves 
to system sizes falling in a region centered around dis- 
crete values Ln ~ 2TTn/qc, with these regions having 
a size of order An. By picking exactly L„ = 2TTn/qc 
(with L < Lthrcsh), so that the characteristic momentum 
mode is exactly resonant with the tunneling rate of zero 
momentum mode atoms, one can maximize the rate of 
MQST decay and maximize the transfer of particles to 
the characteristic modes ±qc before thermalization sets 
in. For the very special case where one picks the smallest 
system size where MQST decay is possible {Li = 2TT/qc) 
we find that the system reaches a steady state only at 
very long times, since only very few system modes par- 
ticipate in the dynamics. 

Interestingly, by taking advantage of the MQST de- 
cay mechanism, the system considered in this paper may 
be used to construct stable states with a sharp momen- 
tum distribution around nonzero characteristic momen- 
tum modes. One simply decouples the two wires after 
enough particles are transfered to the empty wire before 
thermalization sets in, i.e. while the momentum distribu- 
tion is sharp. The decoupled wires are integrable systems 
(see Ref. and the system does not thcrmalizc after 

decoupling. The empty wire in particular has almost all 
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FIG. 13: (Color online) A 3D plot of the momentum distribution in the empty wire vs time and momentum. The parameters 
are shown in the plot and are the same as in Figure [T] One can clearly distinguish three stages in the dynamics corresponding 
to MQST, unstable resonant dynamics producing a sharply peaked momentum distribution, and finally thermalization. 



particles occupying the characteristic modes ±gc- This 
is akin to having large counterpropagating currents in 
the same wire, where interestingly we observe no current 
decay due to the intcgrabiUty of the system. 

We briefly mention here that there exists another way 
of preserving the sharp momentum distribution. Instead 
of decoupling the wires (while the momentum distribu- 
tion is still sharply peaked around characteristic mo- 
menta ±(jfc), one can also achieve the same effect by in- 
troducing a large uniform potential bias V (with V ^ fi) 
between the two wires, the empty wire being at a lower 
potential (if instead the empty is at a higher potential, 
then the self trapping trivially survives and there is not 
much transfer between wires since it is more favorable 
for them to stay in the full wire). The behavior at short 
times is essentially identical to the nonbiased case con- 
sidered in this paper, the only difference being the trivial 
rescaling of the parameters in the biased case. There is 
however a stark difference in behavior between the two 
cases at times greater than 0(1/ J). In this situation, the 
sharp momentum distribution in the empty wire remains 
robust, and MQST survives for long times. The height 



of the peaks in the momentum distribution at first sight 
seems to saturate at the same value as before (i.e. 4.5 
percent of total TV for K/tt = 100 and L = Liq). In 
our numerics, however, we do notice a slow leak of par- 
ticles from the zero momentum mode in the full wire to 
nonzero momentum modes in the empty wire surround- 
ing, and including, the highly populated characteristic 
modes ±qc. This effect slowly increases the number of 
particles in the empty wire, while mostly maintaining 
the sharpness of the momentum distribution. There is 
noticeable broadening developing at later times, and we 
expect that the peaks will once again saturate to a value 
similar to that found in the nonbiased case before ther- 
malization sets in, but this requires further investigation. 
Interestingly we observe essentially no transfer to zero 
momentum modes in the empty wire, in stark contrast 
to what happens with no bias. Because of this slow de- 
cay that persists, we expect that at time scales much 
larger than 1/J (beyond what we have investigated nu- 
merically) MQST eventually is destroyed, but the time 
scales for this seem to be very high. More work needs to 
be done in the biased case. Note once again that in this 
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scenario the wires remain coupled at all times t > 0, the 
only thing we do is introduce an additional potential bias 
before coupling the wires at t = 0. 

As a final remark we point out that in the setup we an- 
alyzed dynamics is particularly interesting for very weak 
but nonzero interactions (see Fig. [T3|) . In this case one 
can clearly see three stages of evolution: (i) macroscopic 
self trapping, (ii) resonant excitation of characteristic 
modes, and (iii) thermalization between the wires. As 
interactions get stronger the second stage gradually dis- 
appears and the sharp momentum distribution at inter- 
mediate times becomes less pronounced. Nevertheless, 
according to our numerical results, all three features can 
be observed for K/tt > 10, which is realistic for present 
experiments. Another important prediction of our anal- 
ysis, which can be easily verified in experiments, is very 
strong and nonmonotonic dependence of the MQST de- 
cay time on the system size. In particular, we find that 
in the interval Lmin < L < Lthresh this decay time is 



relatively short r ~ 1/ J if the resonance condition is sat- 
isfied (i.e. there is an unstable mode) and very long if it 
is not satisfied. This result is valid even in the systems 
with larger interactions, where a narrow momentum dis- 
tribution at intermediate times never forms. 



Acknowledgments 

We acknowledge helpful discussions with R. Barankov 
and A. M. Rey. We also thank I. Bloch for suggest- 
ing this problem and sharing some details on the exper- 
imental system. R.H. acknowledges helpful discussions 
with D. W. Hutchinson on the implementation of TWA 
in continuum systems and acknowledges E. Dalla Torre 
for helpful comments on the manuscript. This work was 
supported by AFOSR YIP and Sloan Foundation. 



[1] S. Giovanazzi, A. Smerzi, and S. Fantoni, Phys. Rev. 
Lett. 84, 4521 (2000). 

[2] S. Levy, E. Lahoud, I. Shomroni, and J. Steinhauer, Na- 
ture 449, 579 (2007). 

[3] I. Zapata, F. Sols, and A. J. Leggett, Phys. Rev. A 57, 
R28 (1998). 

[4] S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy, 

Phys. Rev. A 59, 620 (1999). 
[5] M. Albiez, R. Gati, J. Foiling, S. Hunsmann, M. Gris- 

tiani, and M. K. Oberthaler, Phys. Rev. Lett. 95, 010402 

(2005). 

[6] S. Foiling, S. Trotzky, P. Gheinet, M. Feld, R. Saers, 
A. Widera, T. Miiller, and I. Bloch, Nature 448, 1029 
(2007). 

[7] S. Trotzky, P. Gheinet, S. FoUing, M. Feld, U. Schnor- 

rberger, A. M. Rey, A. Polkovnikov, E. A. Demler, M. D. 

Lukin, and I. Bloch, Science 319, 295 (2008). 
[8] T. Giamarchi, Quantum Physics in One Dimension 

(Glarendon Press, Oxford, 2004). 
[9] A. Widera, S. Trotzky, P. Gheinet, S. Foiling, F. Gerbier, 

I. Bloch, V. Gritsev, M. D. Lukin, and E. Demler, Phys. 

Rev. Lett. 100, 140401 (2008). 
[10] I. Bloch, private communication. 

[11] S. D. Ruber and E. Altman, arXiv:0907.0267 (unpub- 
lished) . 

[12] L. D. Landau and E. M. Lifshitz, Mechanics, Third Edi- 
tion (Pergamon Press, 1976). 

[13] E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963). 

[14] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, 
Phys. Rev. Lett. 98, 050405 (2007). 

[15] M. Rigol, Phys. Rev. Lett. 103, 100403 (2009). 

[16] M. Tabor, Chaos and Integrability in Nonlinear Dynam- 
ics (John Wiley & Sons, 1989). 

[17] V. I. Arnold, Mathematical Methods of Classical Mechan- 
ics (Springer- Ver lag, 1978). 

[18] V. I. Arnold, Russ. Math. Surv. 18, 85 (1963). 

[19] A. N. Kolmogorov, Dokl. Akad. Nauk. SSSR 98, 525 
(1957). 

[20] J. Moser, Nachr. Akad. Wiss. Goettingen Math. Phys. 



Kl, 1 (1962). 

[21] I. E. Mazets, T. Schumm, and J. Schmiedmayer, Phys. 

Rev. Lett. 100, 210403 (2008). 
[22] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 

900 (2006). 
[23] M. Gazalilla, J. Phys. B 37, SI (2004). 
[24] D. Walls and G. Milburn, Quantum Optics (Springer- 

Verlag, Berlin, 1994). 
[25] G. Gardiner and P. ZoUer, Quantum Noise (Springer- 

Verlag, Berlin Heidelberg, 2004), 3rd ed. 
[26] M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, 

S. M. Tan, M. J. GoUett, D. F. Walls, and R.Graham, 

Phys. Rev. A 58, 4824 (1998). 
[27] A. Polkovnikov, Phys. Rev. A 68, 053604 (2003). 
[28] P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh, 

and G. W. Gardiner, Advances in Physics 57, 363 (2008). 
[29] A. Polkovnikov, arXiv:0905.3384 (2009). 
[30] A. Polkovnikov, S. Sachdev, and S. M. Girvin, Phys. Rev. 

A 66, 053607 (2002). 
[31] A. Polkovnikov, Phys. Rev. A 68, 033609 (2003). 
[32] A. Polkovnikov and V. Gritsev, Nature Physics 4, 477 

(2008). 

[33] If two atoms occupying a double well are driven away 
from equilibrium then at small interactions one observes 
single frequency Josephson oscillations while at strong in- 
teractions one observes oscillations described by two fre- 
quencies. The first high frequency small amplitude oscil- 
lation describes single atom virtual tunneling to the un- 
occupied well, and the second, low frequency oscillation, 
describes simultaneous co-tunneling of two atoms. As one 
increases the number of atoms per well the second collec- 
tive tunneling mode becomes exponentially suppressed 
and the first mode describes MQST. 

[34] The precise values of allowed momenta are determined 
by boundary conditions, which we assumed to be pe- 
riodic for the purpose of this paper. Qualitatively the 
picture should not change if we use any other boundary 
conditions. 



